# Packages. 
#xfun::install_github("kuriwaki/ccesMRPprep")
#xfun::install_github("kuriwaki/ccesMRPrun")
xfun::install_github("kuriwaki/ccesMRPviz")
library(ccesMRPprep)
library(ccesMRPrun)
library(tidyverse)
library(ccesMRPviz)
library(ggplot2)
library(dplyr)

# Downloading the CES data. It may take about 30 seconds to download.
ces = get_cces_dataverse("cumulative")
ces20 = get_cces_dataverse("2020")
ces16 = get_cces_dataverse("2016")

# Cleaning the CES data.
ces20_std = ccc_std_demographics(ces %>% filter(year == 2020))
ces16_std = ccc_std_demographics(ces %>% filter(year == 2016))

# Did a candidate or political campaign organization contact you during the [2020/2016] election? <1> Yes <2> No
# There's a follow up (Xb) about the type of 
ces20_question = get_cces_question("CC20_431a", 
                                   dataframe = ces20, 
                                   year = 2020, 
                                   qID = "CC20_431a")
ces16_question = get_cces_question("CC16_425a", 
                                   dataframe = ces16, 
                                   year = 2016, 
                                   qID = "CC16_425a")

# Bringing it all together.
ces20_join  = inner_join(ces20_std %>% mutate(case_id = as.character(case_id)),
                         ces20_question) %>%
  mutate(response = ifelse(response == "Yes",1,
                           ifelse(response == "No",0,NA)))
ces16_join  = inner_join(ces16_std %>% mutate(case_id = as.character(case_id)),
                         ces16_question) %>%
  mutate(response = ifelse(response == "Yes",1,
                           ifelse(response == "No",0,NA)))

fm_brm <- response ~  gender + educ_3 + race + (1|county_fips)

acs_pull_updated <- function (varlist, varlab_df = ccesMRPprep::acscodes_df, year = 2018, 
                              states = NULL, dataset = "acs5", geography = "county") 
{
  acs_df <- get_acs(geography = geography, year = min(max(year, 
                                                          2010), as.numeric(str_sub(Sys.time(), 1, 4))), survey = dataset, 
                    variables = varlist, state = states, geometry = FALSE) %>% 
    filter(!str_detect(.data$NAME, "Puerto Rico")) %>% rename(count = .data$estimate, 
                                                              count_moe = .data$moe) %>% mutate(count = replace_na(.data$count, 
                                                                                                                   0))
  acs_lbl <- select_if(acs_df %>% left_join(varlab_df, by = "variable") %>% 
                         mutate(year = year) %>% 
                         mutate_if(haven::is.labelled, haven::as_factor) %>% 
                         mutate(age = coalesce(.data$age_5, .data$age_10)) %>% 
                         select(acscode = .data$variable, .data$year, .data$GEOID, 
                                matches("(gender|female|age|educ|race)"), .data$count, 
                                .data$count_moe) %>% select(-matches("age_(5|10)")), 
                       ~any(!is.na(.x))) |> 
    rename(county_fips = GEOID)
  acs_lbl
}


acs_tab <- acs_pull_updated(
  varlist = acscodes_sex_educ_race,
  varlab_df = acscodes_df,
  geography = "county",
  dataset = "acs5",
  year = 2019)
acs_tab16 <- acs_pull_updated(
  varlist = acscodes_sex_educ_race,
  varlab_df = acscodes_df,
  geography = "county",
  dataset = "acs5",
  year = 2016)


poststrat <- get_poststrat(
  cleaned_acs = acs_tab, 
  formula = fm_brm) |> 
  filter(county_fips %in% ces20_join$county_fips)
poststrat16 <- get_poststrat(
  cleaned_acs = acs_tab16, 
  formula = fm_brm) |> 
  filter(county_fips %in% ces16_join$county_fips)

ces20_join <- filter(ces20_join, county_fips %in% poststrat$county_fips)
ces216_join <- filter(ces16_join, county_fips %in% poststrat16$county_fips)
# Fitting the model
fit <- fit_brms(fm_brm, ces20_join, verbose = FALSE, .backend = "cmdstanr", .cores = 8)
fit16 <- fit_brms(fm_brm, ces16_join, verbose = FALSE, .backend = "cmdstanr", .cores = 8)

# Generating the estimates. The code assumes respondents in every unit, so it needs to be cleaned a bit
drw <- poststrat_draws(fit, poststrat_tgt = poststrat[poststrat$county_fips
                                                      %in% fit[["data"]]$county_fips, ],
                       area_var = "county_fips")
drw
drw16 <- poststrat_draws(fit16, poststrat_tgt = poststrat16[poststrat16$county_fips
                                                            %in% fit16[["data"]]$county_fips, ],
                         area_var = "county_fips")
drw16

# Looking at the results
mrp_val <- summ_sims(drw, area_var = "county_fips")
mrp_val16 <- summ_sims(drw16, area_var = "county_fips")

saveRDS(mrp_val, "temp/mrp_2020.rds")
saveRDS(mrp_val16, "temp/mrp_2016.rds")
